#!/usr/bin/env python
#coding=utf-8

import sys
sys.path.append("../modules")
import numpy as np
import pylab as P
import numpy as N
import scipy as S
from scipy.signal import detrend
from scipy.stats import mode

import GunnarEEG as eeg_analyser
import eeg_viewer
import os

from __config__ import DATAPATH

s_f = 5.

if __name__ == "__main__":
   
    fname = os.path.join(DATAPATH, 'Gunnar_2010_04_13/Channel_names.csv')
    eeg_analyser.fname = os.path.join(DATAPATH, 
            "Gunnar_2010_04_13/3Gunnar1-export.bin")
    coords_fname = os.path.join(DATAPATH, "Gunnar_2010_04_13/sphere_1005.txt")

    elecnames = eeg_analyser.GetChannelNames(fname)
    RealCoords = eeg_analyser.GetStandardCoordinates(elecnames, coords_fname) 
    channels = N.arange(128,160,1)
    wanted_channels = N.array([0,1,2,3,4,5,6,7,8,9,10
                               ,11,12,13,16,17,18,19,20
                               ,21,22,23,24,25,26,27,28,29,30])
    data = eeg_analyser.ReadBinary(eeg_analyser.fname,
            160,
            channels, 
            readnum_dp='all',
            buffermem=100) * 3.6328e+014
    data = detrend(data, axis=-1)
   
    # find stimuli markers for both hemispheres
    marker = eeg_analyser.GetMarkerPosFromData(data[-1,:])
    marker_diff = mode(N.diff(marker))[0] # find the distance between ipsilateral markers in
    marker = marker[:-1]
    marker2 = marker + (marker_diff/2).astype(int)
    marker_comp = N.sort(N.hstack([marker,marker2]))
    
    data = data[:-1,:]
    
    data = eeg_analyser.InterpolationFromData(data.T, marker_comp, [0,20], mode='dp')
    
    yfilt2 = N.array(
        [eeg_analyser.BPFilter(data[:,k],
            lp=500,
            hp=900,
            sf=5000,
            width=30,
            gpass=5,
            gstop=10) 
         for k in xrange(data.shape[1])]).swapaxes(0,1)
    
    bp2_av1 = eeg_analyser.GetTrialsFromData(yfilt2,marker,[0,1000],mode='dp').mean(axis=2)
    time = np.arange(bp2_av1.shape[0])/s_f

    # calculate scalp map
    timepoint = int(21 * s_f) #timepoint to calc scalp map
    
    InterpCoords_2d, pot1 = eeg_analyser.CalcHemSplineMap(
        RealCoords[wanted_channels,:],
        bp2_av1[timepoint,wanted_channels],
        diameter_samples=30,
        type='2d',
        n=7,
        m=3)

    InterpCoords_2d, csd1 = eeg_analyser.CalcCSDSplineMap(
        RealCoords[wanted_channels,:],
        bp2_av1[timepoint,wanted_channels],
        diameter_samples=30,
        type='2d',
        n=20,
        m=3)

    pot_x,pot_y,pot_z = eeg_analyser.MeshScalpMap(InterpCoords_2d,pot1,gridsize=100)
    csd_x,csd_y,csd_z = eeg_analyser.MeshScalpMap(InterpCoords_2d,csd1,gridsize=100)
    
    P.pcolormesh(pot_x,pot_y,pot_z)
    P.colorbar()

    P.figure()
    P.pcolormesh(csd_x,csd_y,csd_z)
    P.colorbar()
